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Abstract: We study how fluctuations in fluid dynamic fields can be dissipated or amplified within 
the characteristic spatio-temporal structure of a heavy ion collision. The initial conditions for a fluid 
dynamic evolution of heavy ion collisions may contain significant fluctuations in all fluid dynamical 
fields, including the velocity field and its vorticity components. We formulate and analyze the 
theory of local fluctuations around average fluid fields described by Bjorken's model. For conditions 
of laminar flow, when a linearized treatment of the dynamic evolution applies, we discuss explicitly 
how fluctuations of large wave number get dissipated while modes of sufficiently long wave-length 
pass almost unattenuated or can even be amplified. In the opposite case of large Reynold's numbers 
(which is inverse to viscosity), we establish that (after suitable coordinate transformations) the 
dynamics is governed by an evolution equation of non-relativistic Navier-Stokes type that becomes 
essentially two-dimensional at late times. One can then use the theory of Kolmogorov and Kraichnan 
for an explicit characterization of turbulent phenomena in terms of the wave-mode dependence of 
correlations of fluid dynamic fields. We note in particular that fluid dynamic correlations introduce 
characteristic power-law dependences in two-particle correlation functions. 
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1 Introduction 

It is a long-standing idea, first articulated by Landau in the 1950s, that the evolution of matter 
compressed in nuclear collisions lends itself to a fluid dynamical description [1]. On the experimen- 
tal side, characteristic correlations of particle production with the event plane had been interpreted 
as qualitative support for fluid dynamic behavior since the very first relativistic heavy ion colli- 
sion experiments at the BEVALAC in the 1980s [2]. Early qualitative predictions, based on fluid 
dynamics, include notably the argument [3] that the second harmonics of the azimuthal particle dis- 
tribution (elliptic flow f 2) changes at mid-rapidity from out-of-plane to in-plane emission at higher 
center of mass energy. This was confirmed experimentally in heavy ion collisions at the Brookhaven 
National Laboratory (BNL) Alternating Gradient Synchrotron (^/snn < 5 GeV), and at the CERN 
Super Proton Synchrotron ( a /snn < 20 GeV), for a review see e.g. Refs. [4]. Also, fluid dynamic 
arguments provided early on some qualitative understanding of the dependence of elliptic flow on 
particle species, and on the energy and rapidity dependence of the collective sidewards displacement 
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of particle production at projectile rapidity (sidewards flow v\) [4]. However, conclusions remained 
largely limited to the qualitative statement that the observed flow in semi-central collisions "retains 
some signature of the pressure in the high density region created during the initial collision" [5]. 

This changed soon after the start of the Relativistic Heavy Ion Collider (RHIC) in the year 2000, 
when several groups [6-9] noted that fluid dynamic simulations of Au+Au collisions at ^/snn < 200 
GeV can account quantitatively for the main manifestations of collectivity at RHIC, including the 
dominant elliptic flow signal at mid rapidity and its dependencies on transverse momentum, cen- 
trality and particle species. These studies were based on simplified 2+1-dimensional simulations, 
following Bj or ken's argument that the initial conditions for fluid dynamic fields are close to lon- 
gitudinally boost-invariant, and that this boost-invariance is preserved by the fluid dynamics [10]. 
Moreover, early comparisons to RHIC data relied on ideal fluid dynamic equations of motion with- 
out dissipative effects. Soon afterwards, Teaney [11] observed in 2003 that even very small values 
of the ratio rj/s of shear viscosity over entropy density induce dissipative effects that result in a 
sizable reduction of the elliptic flow signal. Therefore, to the extent to which uncertainties in the 
comparison of fluid dynamic simulations with data can be controlled quantitatively, measurements 
of collective flow in heavy ion collisions provide a sensitive tool for constraining transport properties 
of QCD matter. This is one of the main motivations for the development of more and more detailed 
fluid dynamic simulations of relativistic heavy ion collisions in recent years, see e.g. [12-15] for 
recent reviews. 

Ideal fluid dynamics is determined fully by the equation of state and conservation laws. For 
causal viscous fluid dynamics, transport properties and relaxation times enter in addition. But the 
equation of state, transport properties and relaxation times are in principle calculable from first 
principles of a given quantum field theory. Therefore, a comparison of fluid dynamic simulations 
to data of heavy ion collisions has the potential of constraining properties of QCD matter that are 
fundamental in the sense that they are most directly related to the QCD lagrangian. In practice, 
the bottleneck for such a program is the limited control over the initial data that are evolved fluid 
dynamically 1 . Early comparisons of fluid dynamic simulations with RHIC data [6-9] employed a set 
of smooth, event- averaged initial conditions that were specified via the average collision geometry 
in terms of a transverse energy (or entropy) distribution with vanishing flow at initial times. Even 
within this limited set of initial conditions, one observed that differences in the initial transverse 
profile of phenomenologically motivated models could result in variations of the dominant elliptic 
flow signal by up to 30% [16, 17]. 

Within recent years, there has been a growing realization of the importance of event-by-event 
fluctuations in constraining the initial conditions of fluid dynamic evolution in heavy ion colli- 
sions. In particular, event-averaged initial conditions reflect the symmetries of the almond-shaped 
nuclear overlap region of finite impact parameter collisions and can therefore give rise only to 
dipole, quadrupole and higher even moments in the initial density distribution. In marked contrast, 

1 In principle, both the initial conditions for fluid dynamic evolution as well as the conditions for decoupling from 
the fluid dynamic evolution ('freeze-out') imply assumptions that are difficult to control. However, fluid dynamic 
evolution occurs in response to pressure gradients that are much larger at initial times. Therefore, initial conditions 
are thought to be the dominant source of uncertainties. 
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the azimuthal momentum distributions measured by all experiments at the LHC [18-21] and at 
RHIC [22, 23] show a prominent third harmonic moment V3, as well as non-vanishing moments V\ 
and v 5 in addition to the expected even ones. These structures had been attributed previously to 
other speculative effects (" Mach-cone" , "ridge"), but as pointed out first by Alver and Roland [24] 
(see also Sorensen [25] for a related earlier suggestion) they emerge most naturally from the fluid 
dynamical evolution of initial density inhomogeneities. In addition, fluctuations increase the spatial 
eccentricity of initial transverse density distributions, and this accounts naturally for the fact that 
elliptic flow values remain sizable in the most central collisions and for smaller colliding systems [26] . 
There is by now compelling evidence that the dynamical evolution of fluctuating initial conditions 
is a prerequisite for a detailed quantitative understanding of flow in heavy ion collisions [27, 28]. 
And since the various flow moments v n depend differently on the event- averaged initial state and 
its event-by-event fluctuations, analyzing the dynamical evolution of these initial fluctuations pro- 
vides a novel tool for constraining the main uncertainty in fluid dynamical simulations of heavy ion 
collisions. 

There has been a significant effort recently in studying fluctuating initial conditions in heavy 
ion collisions [29-34] and studying their propagation in full fluid dynamic simulations or transport 
models [35-43]. Precursors of these developments include e.g. Refs. [44-46]. The recent efforts 
focussed mainly on initial density inhomogeneities. But more general fluctuating initial conditions 
are conceivable. For instance, (non fluid-dynamical) initial fluctuations in the flow field may 
be expected to accompany a fluctuating initial energy density profile e [37]. Even if fluctuations 
in the initial spatial distribution may turn out to be sufficient to account for the measured flow 
components, it is clearly important to constrain such other conceivable sources of initial fluctuations 
since these may confound any quantitative interpretation of flow phenomena aimed at an extraction 
of rj/s and other fundamental properties of QCD matter. This argues for treating fluctuations in 
all fluid dynamical fields democratically. 

The present paper aims at supplementing the current discussion of fluctuating initial conditions 
with a model study in which the propagation of initial fluctuations can be followed in a very explicit, 
partly analytical way. To this end, we formulate the fluid dynamic evolution of fluctuations in all 
fluid dynamic fields around an event-averaged Bjorken flow profile. The inclusion of fluctuations 
in all fields will provide access to qualitatively novel features such as the dynamical evolution of 
vorticity. It will also allow us to discuss anew how one of the most characteristic manifestations of 
fluid dynamics, namely turbulence, can emerge in the specific expanding geometry of a Bjorken-like 
flow profile. The issue here is not whether heavy ion collisions can display fully developed turbulence: 
It has been pointed out previously (see e.g. Ref. [47]) that the relevant Reynolds numbers are 
typically larger than unity, since they are proportional to the inverse of the normalized viscosity 
rj/s. However, the length and time-scales in heavy ion collisions are so small that Re < (9(100) 
which is well below the conditions under which fully developed turbulence is expected. Rather, 
what is at stake is the suggestion first made by Mishra et al. [48] and further discussed by Mocsy 
and Sorensen [33] that the measurement of the harmonic flow coefficients v n for all values of n 
may provide information about the initial state similar to the power spectrum extracted from 
Cosmic microwave background (CMB) radiation. CMB analysis tools exploit the fact that Hubble 
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expansion dampens vorticity fluctuations, so that the fluid dynamical evolution stays at all time 
scales in a linear non-turbulent regime. It is a priori unclear whether the same situation persists 
for small fluctuations in a Bjorken-type expansion, or whether small fluctuations can become seeds 
of turbulent behavior. Our discussion shall address this question and characterize the limitations 
of a linear treatment of fluid dynamic fluctuations in heavy ion collisions, thus gaining some insight 
into the conditions for onset of turbulent behavior 2 . 

Our work is organized as follows. In section 2, we show that mild extension of models of 
fluctuating initial conditions give rise to fluctuations in all fluid fields. We note in particular that 
fluctuations in the velocity field can have in general a vorticity component as well as a divergent 
component, and that both components may be of comparable size. Furnished with this example 
that fluctuations in all fluid fields may be relevant, we formulate then the equations of motion for 
fluctuations around a Bjorken flow field in section 3, and we solve them in a linearized approximation 
of the evolution in section 4. We then turn in section 5 to the case of turbulent fluctuations, when 
non-linear contributions to the equations of motion matter. In particular, we provide a parametric 
argument that a large class of fluctuating initial conditions around Bjorken flow evolves at late 
times towards an effectively two-dimensional, turbulent system. Motivated by this observation, we 
recall in section 6 pertinent features of turbulence in terms of correlation functions of fluid fields. In 
section 7, we finally relate these remarks to heavy ion phenomenology by showing how correlations 
of fluid dynamic fields enter the one- and identical two-particle correlation functions in a blast wave 
model supplemented with fluctuations. In the conclusion, we finally summarize our main findings 
and provide a short outlook. 

2 Fluctuating initial conditions and vorticity 

We start our discussion with the prototype of an initial density inhomogeneity implemented in 
current event-by-event fluid-dynamical simulations, as used e.g. in Ref. [39]. Fluctuations in the 
initial spatial distribution are described at some initial time r and close to mid-rapidity y = by 
a two-dimensional transverse energy density profile of the form 



Here, the coordinates Xj denote for one specific heavy ion collision the positions of wounded nucleons 
in the transverse plane, as obtained from a Monte Carlo Glauber simulation, see e.g. [53]. A class of 
events corresponds then to a class of independently simulated distributions {xj}, each defining an 
energy density e(x) with different, event-specific fluctuations and each setting the initial data for an 

2 We remark in this context that turbulence may play a role for heavy ion collisions not only in the more narrow 
sense of fluid dynamics as we discuss it in this paper. It has been argued that non-Abelian gauge theories show 
turbulent phenomena such as energy cascades in connection with plasma instabilities and that this could play an 
important role for the time evolution in the early stage of a heavy-ion collision before a hydrodynamical description 
which is based on local thermal equilibrium becomes applicable [49-52]. 




(2.1) 
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individual fluid dynamic evolution. The normalization K in (2.1) can be constrained by data on the 
total transverse energy produced in the collision per unit rapidity [39]. The smearing parameter a 
is a mo del- dependent input that sets the scale of spatial inhomogeneities. Fig. 1 illustrates that this 
model accounts for significant fluctuations in the transverse distribution of wounded nucleons and 
their corresponding energy density (2.1). In Fig. 1, we have chosen a smearing parameter a = 0.4 
fm, consistent with previous simulations of event-by-event fluid dynamics [39]. 
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Figure 1. (Left hand side) The transverse spatial distribution of nucleons as obtained from a Monte Carlo 
Glauber simulation of a Pb-Pb collision at 6 fm impact parameter. Nucleons of the two colliding nuclei 
are characterized by red and black circles, respectively. The radii correspond to a black sphere inelastic 
nucleon-nucleon cross section of 63 mb at snn = 2.76 TeV. 'Wounded' nucleons that interact are denoted 
by smaller blue and green circles, respectively. (Right hand side) Initial distribution of transverse energy 
density (2.1) corresponding to the distribution of wounded nucleons on the left hand side. This plot is in 
arbitrary units. 



One possibility is to initialize fluid dynamical fields at initial time r with event- wise fluctuations 
in the transverse energy density e(x), but with an exactly vanishing non-fluctuating flow field in 
the two transverse directions 1 and 2 at initial time To, 

^(x^tq) = w 2 (x,y,Tb) = 0, (2.2) 

and in the rapidity component of the flow vector 

u»(x,j/,t„) = 0. (2.3) 

However, a larger class of initial conditions is conceivable, since the initial state may also display 
fluctuations in the initial flow fields, and since both energy density (2.1) and the flow fields could 
depend on rapidity y. Fluctuations in the initial velocity fields have been discussed recently e.g. 
within a model of the early (Non-Equilibrium) dynamics based on free streaming [38]. For more 
discussions of initial flow and its influence on HBT radii see also Refs. [54-58]. 

Some fluctuations in the flow field will be generated by the fluid dynamic response to initial 
density fluctuations and may therefore be regarded as being implicitly included in the ansatz (2.1) 
- (2.3). However, such fluid dynamically generated flow fluctuations will be constraint in scale 
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and size to the fluctuations in energy density, and they will lack by construction some qualitative 
features of the most general fluctuating flow field, such as vorticity. Vorticity 3 characterizes the 
solenoidal part of a general three-dimensional flow field u^(x), 

I \ (d 2 u y - d y u 2 )\ 
Uj = (Curl u)j = \ {dyUi — diUy) . (2.4) 
V <9i« 2 - d 2 ui J 

We note that the fluid dynamical evolution of fluctuations in vorticity and energy density decouples 
as long as these fluctuations are small and can be treated in a linearized evolution (see section 4). 
This is in marked contrast to fluctuations in the irrotational part of u 1 (sound modes) that can be 
driven by fluctuations in energy density. This indicates that one cannot expect to generate sizable 
values of vorticity by evolving initial conditions of the form (2.1)- (2.3). However, if fluctuations in 
vorticity are part of the initial conditions, then they will propagate and may display particularly 
interesting dynamical features, as discussed in sections 4 and 6, respectively. 

To set the stage for our discussion in later section, we demonstrate now that relatively mild 
extensions of (2.1) can lead naturally to fluctuations in velocity, including a non- vanishing solenoidal 
component (2.4). It is one arguably mild extension of (2.1) to associate the transverse region around 
a single wounded nucleon not with an energy density, but with an energy- momentum tensor T£ u , 
such that the initial energy-momentum tensor of the entire nucleus-nucleus collision takes the form 

iVpart 

T^(T,x\x 2 ,y)= ^THr^-oix 2 -^,?/). (2.5) 
i=i 

In general, equation (2.5) can account for non-vanishing fluctuating initial conditions in both en- 
ergy density e(x) and flow u^(x). For instance, neglecting for simplicity non-ideal, shear viscous 
contributions to (2.5) and assuming an ideal equation of state e(x) = 3p(x), one can write the 
initial conditions for e and u 1 at some fixed rapidity y and initial time tq to linear order in -u- 7 in 
the form 

e(x\x 2 ) (l, |u*(s)) = T^{r Q ,x\x 2 ,y). (2.6) 

The transverse energy density associated to a single wounded nucleon is given by e w (x) = T°°(x). 
Therefore, the 0-component of (2.6) defines an equation of the type (2.1) for the energy density, 

ATpart 

e(x\x 2 ) = ^e^-x^-x 2 ), (2.7) 

i=i 

but it has also a fluctuating initial flow field defined by the spatial components of (2.6), 



V^part , I _ I 2 _ „2\ j ( I _ 1 2 _ 2\ 

T^ru^-xj^x^-x 2 ) 



Here, we deviate from the standard notation of vorticity in terms of a cartesian three-vector by adopting the 
notation of vorticity to light cone coordinates, u J = (it 1 ,u 2 ,u y ). We note that the three components of vorticity 
(ll>i,u>2, W3) do not form the spatial part of a four- vector and it makes no sense to contract them the spatial part the 
metric. 
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Figure 2. Velocity distribution corresponding to Fig. 1 if one assumes for the contribution of each wounded 
nucleon a random velocity in the transverse plane drawn from a Gaussian distribution of width (|v|) = 0.1c. 
Shown is the innermost area —3.5 fm < x 1 , x 2 < 3.5 fm of the transverse plane. 



Here u 3 w is defined by writing equation (2.6) for the energy momentum tensor associated to a single 
wounded nucleon. The size of initial fluctuations in the velocity field (2.8) depends on how the 
initial transverse motion associated to the the single wounded nucleon is modeled. Taking guidance 
from blast wave models, one may choose for u J w e.g. an azimuthally symmetric radial flow field 

j j 

with some radial dependence w, «i,( x ) = |^-x-'| w (l x ~~ x «l) f° r 3 ' — 1j2 and u y = 0, say. For such 
an ansatz, one checks easily that (2.8) defines in general a flow field of non-vanishing transverse 
divergence and non- vanishing longitudinal vorticity, 

cW(x) + «9 2 u 2 (x)^0, 
W3 (x)^0. 

A more general ansatz may be based on the observation that in general, the transverse energy 
deposited by a single wounded nucleon in a finite window of rapidity recoils against transverse 
momentum outside this rapidity window. This argues for a net transverse velocity component Vj 
associated to the contribution of each wounded nucleon in (2.5). To illustrate this effect, we assume 
that each wounded nucleon in the sample shown in Fig. 1 is associated with a non-relativisticaliy 
small random transverse velocity component Vj, drawn from a Gaussian distribution of width (|v|) = 
0.1c. The resulting initial transverse flow field (2.8) is shown in Fig. 2. By comparing to Fig. 14 
of Ref. [38], we note that a seemingly comparable flow field can be obtained in a model of early 
non-equilibrium dynamics based on free streaming where the contribution of every particle is taken 
to be deiocalized in position space over a small volume. 
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Figure 3. The left hand side shows the absolute value of vorticity \d\u 2 — d2U l \ for the velocity field shown 
in Fig. 2. Similarly, the right hand side shows the absolute value of the divergence of the fluid velocity 
\d\u l + d2U 2 \ for the same velocity field. The color coding is the same on both sides. 

In Fig. 3, we plot the absolute value of the longitudinal vorticity \u)?\ = \diu 2 — c^w 1 ] and the 
transverse divergence \diii 1 + c?2M 2 | for the flow field shown in Fig. 2. Inspection of this figure shows 
that both components fluctuate with a similar magnitude and over similar transverse length scales. 

From the simple models discussed in this section, we conclude that fluctuating initial conditions 
may comprise in general fluctuations in both the energy density and the velocity fields, and that 
fluctuations in the irrotational and solenoidal components of the velocity fields may be of comparable 
size. Lacking a general argument for the relative size of fluctuations in m j and e will motivate us to 
treat all conceivable sources of event-by-event fluctuations on a common footing in the subsequent 
discussion. 

3 Fluid dynamic equations of motion for relativist ic heavy ion collisions 

On time and length scales that are large compared to the typical relaxation times and lengths for 
thermal and chemical equilibration, relativistic fluid dynamics provides an effective theory for the 
multi-particle system produced in heavy ion collisions. A large set of experimental observations 
support the assumption that in ultra-relativistic heavy ion collisions the range of validity of this 
effective fluid dynamical description is large and comprises bulk hadron production up to a few GeV 
in transverse momentum [12-14]. In this section, we recall first shortly the fluid dynamic equations 
of motion. We focus then on the Bjorken model that defines a particularly simple expanding 
geometry and that encodes important features of a relativistic heavy ion collision. Regarding the 
Bjorken model as defining the average fluid dynamical field, we then discuss how fluctuations in 
energy density and flow propagate in the expanding geometry of a relativistic heavy ion collision. 
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The relativistic hydro dynamic equations for a fluid without any conserved charges read [12-14] 

Dt + (e + p)8X + ^ w A*d a u v = 0, 
(e + p)Du Q + A^dpp - A a v d^ v = 0. 

Here e is the energy density and p is the pressure in the fluid rest frame, tc^ is the viscous part 
of the energy- momentum tensor in the Landau frame, u fl 7i fiu = 0. The partial derivative d a must 
be replaced by the covariant derivative V a if one works with coordinates other than cartesian. We 
work with a cartesian metric of signature g^ u = diag(— 1, 1, 1, 1). The matrix A^ v projects to the 
subspace orthogonal to the fluid velocity, A^ v = g^ u + u^u v . The derivative in the direction of the 
fluid motion is D = w M <9 M . 

The viscous part of the energy-momentum tensor can be expanded in a derivative expansion. 
To lowest order it vanishes, leading to ideal hydrodynamics. The first order contains shear and bulk 
viscocity terms 

tt^ = -277(7^ + (A^V a u a } (3.2) 

where 

o» v = -{A\V a u u + A\W) - -A^(V a u a ) (3.3) 

is transverse (orthogonal to w M ) and traceless. 

To second order in the gradient expansion, the fluid dynamic equations of motion contain 
various relaxation time corrections. It is a peculiar feature of a second order approximation that 
the evolution equations are hyperbolic and propagation is limited to the forward light cone even for 
perturbations of large wave-vector k. For this reason, second order fluid dynamics is often referred 
to as causal viscous fluid dynamics. However, this wanted feature of causality is not guaranteed 
to persist in higher orders of the gradient expansion. More generally, fluid dynamics is a long 
distance effective theory that by its very construction cannot be expected to be reliable for large 
wave- vectors. For the propagation of fluctuations with small gradients (i.e. small wave- vectors k), 
second order fluid dynamics will make only small corrections to a first order treatment. For this 
reasons and to keep the formalism simple, we restrict the discussion in the present paper to first 
order fluid dynamics. For the case of vanishing bulk viscocity (, the corresponding equations of 
motion read 

Z)6 + (e + p)V/-2^, i y i/ = 0. (3-4) 
(e + p)Du v + A^V M p - 2 v A v a V^ a = . (3.5) 

3.1 The Bjorken model 

We are interested in studying fluid dynamic fluctuations in the expanding geometry of a relativistic 
heavy ion collision. The Bjorken model is arguably the simplest formulation of a corresponding 
expanding geometry. Motivated by the idea that in nuclear collisions at ultra-relativistic energy, 
particle production is almost flat in rapidity, Bjorken [10] proposed to formulate initial conditions 
for fluid dynamic fields that are independent of space-time rapidity y = arctanh(x3/|xo|), that 
means e(r, x, y) = e(r, x) and = (a/1 + u 2 , u, u y ) = (vT+ u 2 , u, 0) . If this condition is satisfied 
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(3.1) 



at some initial proper time r = tq, then it persists for all proper times r = a/xq — x\ throughout 
the evolution. This renders the longitudinal evolution trivial, and the numerical task simplifies to 
the solution of a (2+l)-dimensional problem. 

Under the further simplifying assumption that the initial transverse flow field vanishes at initial 
times and that transverse gradients in energy density are absent, the Bjorken model reduces to an 
effectively (l+l)-dimensional toy model that allows for an explicit analytical treatment. In this 
case, the evolution equation for the energy density becomes 



d T e+^- V ^ = 0, (3.6) 



where energy density and pressure are related by the equation of state e = e(p). For what follows, 
it will be useful to rewrite this equation in terms of the enthalpy 

w = e + p = sT , (3.7) 



v — r)/w . (3.8) 



d T e + ™(l- A /)=0. (3.9) 



and the kinematic viscosity 
Eq. (3.6) reads then 

d„e. + . 

t \ 3r / 

For large times r > 1/T and for a small normalized shear viscosity t]/s <C 1, one can replace the 
bracket in (3.9) by unity. We shall be able to employ this approximation in the following. For an 
ideal equation of state e = 3p, one then finds the characteristic time dependencies of the Bjorken 
model for energy density 

/Tn\ 4 /3 

eB j (r)=e Bj (ro)(^) , (3.10) 

and temperature 

TB j (r)=T Bj (r )f^ N ) 1/3 . (3.11) 



T 

For a time- independent normalized viscosity rj/s, the ratio 

^Bj(r) = ^Bj(ro) /Tq\2/3 

T T V T J 

decreases. Therefore, replacing the bracket in Eq. (3.9) by unity is an approximation that is 
consistent with the late time behavior. 

We note at this point that in situations with strong (non-Gaussian) fluctuations, the evolution 
equation for averaged fields such as energy density e(r) = (e(r)) gets modified by additional terms, 
see the discussion in Sect. 6.1. For the present paper we assume that these modifications are small 
and can be neglected. 
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3.2 Fluctuations on top of a Bjorken background field 



In this section, we formulate the theory of the dynamics of fluctuations on top of a Bjorken back- 
ground field without transverse gradients. That means that the hydrodynamical fields u^, e when 
averaged over many events follow a Bjorken type solution. However, locally and for a particular 
event we expect deviations which we want to investigate in more detail. We have chosen a Bjorken 
background field for our study mainly for two reasons. First, the analytical simplicity of this back- 
ground will allow for a particularly explicit discussion. Second, the Bjorken model contains essential 
features of realistic expansion scenarios of relativistic heavy ion collisions. 4 

We denote fluctuations on top of the Bjorken flow u^. = (1, 0, 0, 0) by relaxing the constraints 
(2.2) and (2.3) and allowing for local fluctuations in the transverse and rapidity components, 
it 1 , it 2 , vfl . The normalization condition u^u^ = —1 of the local fluid velocity = (u T , w 1 ,u 2 ,u y ) 
implies then 



Here and in what follows, we work in light-cone coordinates r, x±,X2, y with metric g^ u = diag(— 1, 1, 
1, 1/t 2 ). The latin index j is summed over 1,2, y and the corresponding three-dimensional metric 
reads g iJ = diag(l, 1, 1/t 2 ). We consider small local fluctuations in the sense that UjU^x) <C 1. 

Neglecting terms that are parametrically suppressed due to UjU^x) C 1 or due to vjr <C 1, 
we find from Eq. (3.4) and (3.5) the following equations governing the velocities in the transverse 
plane (j = 1,2) and in rapidity direction [j = y) 



Here, the first two terms describe the change in the velocity along the direction of the fluid motion. 
They result from writing Duj = u^d^Uj for small deviations from the Bjorken background. The 
terms in the first square bracket account for two effects. One is the acceleration of the fluid due to 
pressure gradients in the transverse direction. The second term proportional to Uj is dominated by 
the decrease of pressure for increasing r. This dilution of the fluid leads to an acceleration in the 
direction of Uj. Finally, there are effects of viscosity that are similar to the corresponding term in 
the (non-relativistic) Navier-Stokes equation. 

4 It has been pointed out repeatedly that despite its simplicity, the Bjorken model without transverse gradients 
retains important features of the early time dynamics of heavy ion collisions. The argument is based on the observation 
that event-averaged initial energy density distributions show typically only small transverse gradients in the central 
region of the transverse plane; the central region may indeed be approximated by the ansatz u = 0, e(x) = e. 
The transverse evolution of this initial condition may then be thought of qualitatively as being dominated by a 
rarefaction wave that moves from the outside (vacuum) to more and more central positions in the transverse plane at 
late times. At a given position in the transverse plane, the dynamical evolution may be viewed as being characterized 
by the effectively (1+1)- dimensional Bjorken model up to the later time at which the rarefaction wave reaches the 
corresponding transverse position. Based on such considerations, the Bjorken estimates for the time-dependence of 
energy-density (3.10) and temperature (3.11) are used regularily in simple phenomenological estimates. 




(3.13) 



d T Uj + u l diUj H — [djp + Uj(d T p + u l dip)\ — v -djdiU 1 + did l Uj = . 



(3.14) 
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In addition to eq.(3.14), one finds under the same assumptions for small local fluctuations 
around the Bjorken background the equation of motion for the internal energy density 



d T e + u 3 dj€ + w 







- + d jU j 


-v 


T 



diUjd l u j + diUjd j u l — -ditfdjV? 



0. 



(3.15) 



Here, the first two terms describe the change along the fluid direction of motion. The first square 
bracket describes dilution effects; the first term ~ - is due to the expansion of the Bjorken- 
background in the longitudinal direction while the second term measures the effect of a possible 
dilution (or compression) in the transverse and rapidity directions. Viscous correction that are 
parametrically suppressed due to t)/(wt) <C 1 have been dropped, and the remaining dissipative 
contribution to the evolution of internal energy are given in the last bracket of (3.15). They describe 
how kinetic energy is transferred from the macroscopic motion of the fluid to internal energy. 

It will turn out to be useful to rewrite Eqs. (3.14) and (3.15) in terms of rescaled fluctuations 
in velocity, 

/ T 



and for a rescaled time variable 



3r 4 / 3 



4r, 



1/3 : 



T 



dt 



1/3 



J 3 ' 



(3.16) 



(3.17) 



(Of course, this rescaled time t is not the time variable x° in the laboratory frame.) In what follows, 
we also absorb deviations from the r-dependence of Bjorken's energy density (3.10) in terms of the 
quantity 

-) d, d = \n — — . 

Finally, we assume that the shear viscosity v follows the Bjorken behavior (3.12). That means, we 
neglect local fluctuations in the kinematic viscosity since they are expected to have only a minor 
effect. The kinematic viscosity z/q can then be written as 



cL 



(3.18) 



v 



1/3 



^Bj(ro). 



(3.19) 



For an ideal equation of state e = 3p, using ^dp = -^^dT = j;dT, and neglecting a dissipation 
term ~ v that is of higher power in the velocity field, one can show that Eqs. (3.14) - (3.15) lead 
to the equation for the (rescaled) velocity (j = 1,2,?/) 



2 l 



(3.20) 
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and for the quantity d 7 



1 ^ 11 2/3 

d t d T + —d T + ^2 v m d m d T + —Vyd y d T + - i? - ^ 



m=l 



2j (^mfn + d n V m )(d m V n + <9„ 



m,n=l 



2 x - 4 

+ ~^ 2 + d mVy)(dyV m + d m Vy) + —(9; 



,2 _ ^2 



m=l 



0. 



(3.21) 



Here and in what follows, we denote the expansion scalar dj of the rescaled velocity fields by 



1 



■d = dxvt + d 2 v 2 + -^d y Vy. 



(3.22) 



In the following we will use both the representation of the equations of motion in Eq. (3.14), (3.15) 
and the one in (3.20), (3.21). In particular the discussion of linear fluctuations in section 4 will be 
largely based on (3.14) and (3.15), while for the discussion of non-linear fluctuations in section 5 
the representation (3.20), (3.21) will be more appropriate. 



4 Linear fluctuations 



In this section we discussion the evolution of fluid dynamical fluctuations that are small enough 
to neglect non-linear terms in (3.14) and (3.15). In addition, it is assumed that the deviation of 
the temperature field from the homogeneous background is small, d < 1. The resulting linearized 
equations describe laminar flow. They apply to systems with sufficiently small Reynolds number, 
as we shall discuss in section 5. 

For a fixed time r and given spatial boundary conditions, one can divide the velocity field 
uniquely into a solenoidal part with vanishing divergence, and an irrotational part with vanishing 
curl, 

Uj = Uj +Uj , 

Divu s = 0, (4.1) 
Curl u 1 = , 

where Div u = djU^ and Curl u is defined in Eq. (2.4). We recall that the derivative operators 
in (4.1) introduce an explicit r-dependence since they involve the three-dimensional metric = 
diag(l, 1, t 2 ). Therefore, the splitting of Uj into a solenoidal and an irrotational part does not 
commute with the r-derivative. 

The field can be represented by the expansion scalar (3.22) 
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It will be convenient to write this expansion scalar as a sum of a transverse and a longitudinal 
contribution, 9 = 9t + 9 y . 

From (3.14), (3.15) we obtain the linearized equations 

d T 6 T - ^0t + (8? + d 2 2 )d - v [\{dl + d 2 2 )6 + {dl + d 2 2 + ±d 2 y )9 T ] = 0, (4.3) 

dr9 y + ^9 y + - i/ [\^d 2 y 6 + (d 2 + dl + ±d 2 )9 y ] = 0, (4.4) 

d T d+\e = 0, (4.5) 

- ^ - + d 2 2 + i^ 2 )^- = 0. (4.6) 

Here, we recall that the quantity d denotes the logarithmic temperature (3.18). The symbol hj in 
(4.6) takes the values hi = /12 = — 2 and /13 = 1. 

Interestingly, the vorticity modes Uj decouple from the velocity divergence 9, the logarithmic 
temperature field d and from each other. Using Fourier decomposition with respect to the spatial 
argument, 

/d 3 k 
— c^(r, k x , k 2 , k y ) e i(fcx»i+fa^+fc»») , (4.7) 

their diffusion-type equation of motion can be directly solved 

To 



u^^A) = Wi (7 b ,Jfe 1> ifc a ,ife3) x ( — V"' e-^^^^U-^. (4.8) 



We assume here a constant z/ as defined in (3.19) and we use (3.17). One sees from equation 
(4.8) that vorticity modes for essentially all wave vectors are dominated at late times by an expo- 
nentially decaying function with a decay time set by the product of kinematic viscosity and the 
square of the wave vector. A somewhat unusal case is the time evolution of modes with k y ^ 
where the exponential damping term is modified by a term oc 1/ In particular, for k\ + k\ = 
but k y ^ the vorticities do not decay exponentially for r — > 00. In addition, the exponential 
decay is modified by a term that decreases algebraically for the transverse components u± and u>2 
and that increases in the longitudinal components U3. For finite times, the algebraic increase of 0J3 
can overcome the exponential decay with viscosity. 

In Fig. 4, we plot the solution (4.8) of the linearized fluid dynamic equations of motion for 
phenomenologically motivated input values, namely a small normalized shear viscosity rj/s = 1/4tt 
and a temperature of 500 MeV at initial time tq = 1 fm/c. This translates into an initial kinematic 
viscosity u ~ 0.03 fm. Most generally, Fig. 4 illustrates the interplay between an exponential decay 
set by kinematic viscosity, and the characteristic algebraic dependencies of the transverse and 
longitudinal vorticity components. More specifically, we have chosen in Fig. 4 wave vectors that 
correspond to fluctuations on length scales between 1 fm and 0.25 fm, as may be regarded as realistic 
for the initial state of the system created in heavy ion collisions. We observe from Fig. 4 that such 
fluctuations are modified but persist over time scales of O(10 fm/c) typical for the expansion history 



- 14 - 



transverse propagation (k2=ky=0) 



s 1.5 





1 1 


Cflj forkj = 1/fm 

CO3 fork j = 2/fm 




. _ . _ . CO3 for kj = 4/fm 




Cflj forkj = 1/fm 




. _ . . _ coj for kj = 2/fm 


- 


_ . cflj forkj = 4/fm 






- 


/ v. 
\\ 




- K K 


s. — 






, 1 , i- -,~;ir--, ■■ 



8 10 



longitudinal propagation (kj=k2=0) 



r<Vj fnr V =1 


1 1 1 


CO3 for ky = 2 




. _ . _ . CO3 for ky = 4 




coj forky = 1 




. _ . . _ C0j forky = 2 


_ / ~ 


_ . coj forky = 4 




s 




y 














/ ^ — 








, 1,1,1 


, 1 , 



10 



T [fm] 



Figure 4. Time dependence of the normalized transverse (i = 1) and longitudinal (i = 3) vorticity 
amplitude (4.8) for modes with wave vectors k\, k<i = k y = in the transverse direction (left hand side) 
and with wave vectors k y , k\ = hi = in the longitudinal direction (right hand side). Input values are 
T(r = lfm/c) = 500 MeV and r]/s = l/4vr. 



of relativistic heavy ion collisions. The figure illustrates that over times scales relevant for the fluid 
dynamic expansion of heavy ion collisions, some fluctuations on phenomenologically relevant scales 
may get amplified rather than dampened. Moreover, the relative attenuations (or amplifications) 
of vorticity components over times of O(10 fm/c) are - within the phenomenologically relevant 
parameter range - very sensitive to the length scale 1 / kj of the fluctuations. From inspection of Eq. 
(4.8), it is also evident that there is a similar sensitivity to the precise choice of the viscosity. 

In addition to the evolution equations for vorticity, there are equations for 9t, y and d that 
we discuss now. These describe sound waves and are best solved in Fourier space. We concentrate 
first on a wave traveling in the transverse direction xi, corresponding to k\ 7^ 0, k% = k y = 0. In 
this case, eq.(4.4) decouples from the others 

d T 9 y + ^-9 y + vk\Q y = (4.9) 
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and can be integrated immediately, 

^(^,0,0) = 9 y (T ,h,0,0) ^) 5/3 e -^f(^o), (4.10) 
Equations (4.3) and (4.5) are coupled, 

d T e T - i- r e T - k\d + vk\ (p T + \e y ) = o, n 

d T d+l(6 T + 6 y ) = 0, 

and depend also on the solution for 9 y . Concentrating for simplicity on the case 9 y = and 
eliminating d, one finds the second order differential equation 

drOr + (~ + \vk 2 ^ d T 6 T + f-L + ±lA 6t = 0. (4.12) 

(We have dropped a term suppressed due to vjr <C 1 in the second bracket.) For vanishing viscosity, 
this equation can be solved in terms of Bessel functions. The two linear independent solutions are 

tV M%) and -MtO- (4 - i3) 

One finds an oscillating behavior with the amplitude increasing algebraically with time proportional 
to (r/r ) 1/6 . For non-vanishing viscosity v there is also an exponential decay which is larger for 
large wavevectors k\. For small values of k\T the two independent solutions of (4.12) are 

proportional to r and r 1 / 3 , respectively. Eq. (4.11) implies that the temperature field d grows 
according to r 2 and r 4 / 3 for these two solutions. For k\ ^ and late enough times r, the solutions 
in (4.13) always have an oscillating behavior, however. For large wave vector k\ one can neglect the 
terms ~ - and ~ \ in (4.12). Up to viscous damping, the solution corresponds to a perturbation 
that propagates with the velocity of sound cs = a/1/3 into the transverse direction. 

We show the time evolution governed by (4.12) on the left hand side of Fig. 5. In close similarity 
to the case of vorticity, we observe that also fluctuations in the transverse velocity divergence 
can persist over time scales relevant in heavy ion collisions. For sound waves in the transverse 
direction, the fluid acts like an efficient low-pass filter, allowing for the unattenuated (or even 
slightly amplified) passage of fluctuations on sufficiently large scales l/k\ > 1 fm, while filtering 
out fluctuations on smaller length scales 1/ki < 0.25 fm. 

We finally turn to sound waves traveling in the rapidity direction y, i. e. k y ^ 0, ki = k 2 = 0. 
Now, equation (4.3) decouples 

dr e T -±8 T + v±kl6 T = (4.14) 

and can be integrated immediately, 

9 T (r, 0, 0, ky) = 9 t (t , 0, 0, k y ) (^) '* e ^4o*Hv^vt) . (4.15) 
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Figure 5. (Left hand side) Velocity divergence amplitude 9t(t)/6t(^^-/c) for sound waves of various 
wave vectors ki, traveling in the transverse direction x\. The r-dependence is calculated from (4.12) for 
T(lfm/c) = 500 MeV and rj/s = 1/47T. (Right hand side) Velocity divergence amplitude 9 y (r)/9y(lim/c) 
for a sound wave traveling in the rapidity direction y, calculated from (4.17) for the same input values. 



The equations for 9 y and d are coupled and depend on 9 T , 

d T e y + ±e v - ±k 2 y d + v±k 2 y (le y + \e T ) = , 

d T d + l(e T + e y ) = o. 

Concentrating on 9? = 0, this yields the following second order differential equation for 9 y (we drop 
again a term suppressed due to vjr <C 1), 

d% + (A + v ±hj) d T 9 y + (- A + ^ 9y = 0. (4.17) 

The two linear independent solutions to this equation for v = are 
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For ky <C 16/3 this becomes r and r~ 5 / 3 . Equation (4.16) implies that this corresponds to perturba- 
tions in the temperature field d that grow like r 2 or decay like r~ 2//3 , respectively. For k 2 > 16/3 the 
solutions (4.17) correspond to an oscillation with a period that increases as a function of time and 
an additional decrease in the amplitude. On the right hand side of Fig. 5, we plot this behaviour 
of 8 y for small non-vanishing viscosity, when the algebraic increase is modified by an exponential 
decay. We observe again that fluctuations of modes with wave vectors k y ~ 0(1) can persist or 
can be amplified over time scales commensurate with the expected expansion duration of heavy ion 
collisions. 

In the limit of large ky, one can translate (4.18) back to position space and one finds that it 
corresponds to an excitation that propagates in the rapidity direction according to 

In general, for a sound mode propagating fastly into an arbitrary direction described by a large 
wavevector k = (k\, fc 2 , k y ), we expect the propagation velocity (^r, ^p, §jr) to satisfy 

This condition is satisfied by the sound modes in 9 y and 9t discussed here. 



5 Turbulent fluctuations 

In full generality, the non-linear equations (3.14) and (3.15) or, equivalently, (3.20) and (3.21) are 
difficult to analyze. One can always use a splitting of the velocity into a solenoidal part that carries 
vorticity u>j and an irrotational part described by the divergence 9. The nonlinear terms in the 
equation of motion will lead to couplings between these fields and to the logarithmic temperature 
field d, however. Intuitively, one expects that perturbations in the fluid divergence propagate quickly 
in the medium with the characteristic velocity given by the velocity of sound. The fluid velocity 
(u l ,u 2 ,u y ) can be small compared to sound propagation. This is often characterized by a small 
Mach number 

Ms= veff™ <<l (51) 

cs 

For the description of the part of the fluid velocity that carries vorticity one can often assume in 
this case a vanishing divergence, 9 = 0, since the fast sound modes can be viewed as decoupled from 
the slower modes dominating the solenoidal part of the fluid velocity. In the present section we will 
study the equations of motion in this situation and show that there are some interesting parallels 
to non-relativistic, incompressible fluids. 

For our discussion in this section, it will be useful to work with the rescaled velocities v,i 
introduced in section 3. For 9 = i? = 0, Eq. (3.20) simplifies then to 



dtV J + VmdmV i + \ v v d y v j + d i d ~ u o (df + d\ + \d 2 y j vj = 0. 

m=l T V T / 



(5.2) 
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Due to the the solenoidal constraint 

d\Vi + d 2 v 2 + ^dyVy = (5.3) 

the temperature field d is not independent of the velocity field. More specific, by taking the 
divergence of (5.2) one derives 

d l + d l + ^ d l) d + 2 (dmv n )(d n v m ) + ^ 2(a m %)(^m) + ^(<9 y ^) 2 = 0. (5.4) 

' m,n=l m=l 

This is an instant of the Poisson equation for d. For given boundary conditions it can be inverted 
to yield d as a (non-local) functional of the velocity field. 

Equation (5.2) has some interesting features 5 . It takes the form of a two-dimensional Navier- 
Stokes equation in situations where v y = and where the dependence of vi, v 2 on y can be neglected. 
Moreover, for a large class of initial conditions at time r = To, the evolution becomes effectively 
two-dimensional for late times t/tq ^> 1. Indeed, both the non- linear velocity term that couples v y 
to Vi and v 2 and the damping term involving the derivative with respect to rapidity contain factors 
that decrease as 1/t 2 . Similarly, the solenoidal constraint (5.3) assumes its two-dimensional form 
in that limit. 

Motivated by the observation that for the particular expansion geometry of the Bjorken model, 
fluctuations are governed by an evolution equation that reduces at late times to a non-relativistic 
Navier-Stokes equation (in rescaled time coordinates), we now discuss the conditions for a non-linear 
turbulent evolution of fluctuations by the very concepts that have proven useful in the classification 
of solutions of the Navier-Stokes equation. To this end, we consider situations where the velocities 
change notably on distances of order I in the transverse direction or for rapidity differences Ay. 
The damping term in Eq. (5.2) leads then to a damping rate (inverse relaxation time) that can be 
characterized in terms of the dimensionless number 

I 2 

(5.6) 



r 2 Ay 2 

This damping rate is of order vq/1 2 for k < 1, and it is of order Vo/{r 2 Ay 2 ) for k ^> 1. For 
characteristic velocities vt in transverse, respectively v y in rapidity direction, the flow can then be 



5 We mention as an aside that in addition to the obvious rotation symmetry in transverse direction and the 
translational symmetries in transverse and rapidity directions, Eqs. (5.2) and (5.3) have also the following space-time 
symmetry 

x m -> x m + V m t (m = 1, 2), 
V^V-2V V ±, (5.5) 



Vi — > Vi 



for constant velocity Vj. Indeed, contributions from the time derivative and from the non- linear advection terms 
cancel. For V y = this corresponds to Galilean symmetry in the transverse plane while the situation is more 
complicated for V y 7^ 0. However, it is not so clear whether this invariance is very useful since there is no translational 
symmetry with respect to time. 
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characterized in terms of the Reynolds numbers 

^, HeW = 
v Q ' vq Ay t 2 



RePl = ^, Re^) = ^l for «<1 (5.7) 



and 

Re (T) = ^TW Re (,) = !^ for R>1 . (5 . 8) 
I/ / ^0 

Obviously, these definitions can be extended to intermediate k, as well. 

If both Reynolds numbers are small, Re^ <C 1, Re^ <C 1, the resulting flow pattern is 
expected to be laminar. The viscous damping term dominates then over the nonlinear terms in Eq. 
(5.2) and velocities are expected to follow a regular behavior dominated by an exponential decay 
in time. More specifically, if the non-linear term in Eq. (5.2) can be neglected one can use Fourier 
decomposition with respect to the spatial arguments and one finds the solution 



h, k 2 , ky) = h, k 2 , h)e ~^ + ^ + ^<Vr A) . 



(5.9) 



This resembles closely the behavior of vorticity in Eq. (4.8). 

In contrast, if both Reynolds numbers are large, Re (T) > 1, Re {y) > 1, one expects a turbulent 
regime. The nonlinear terms dominate now over the viscous damping term in Eq. (5.2) and the 
velocities will change rather irregularly from point to point both in the transverse plane and for 
different values of the rapidity variable y. 

We consider next the case Re^ <C 1 and Re^ ^> 1. For k ^> 1 the derivatives with respect to 
X\ and X2 effectively drop out from Eq. (5.2). One might expect a sort of one-dimensional turbulent 
behavior. However, the unusual explicit time dependence of the non-linear and damping terms 
could spoil this conclusion. Also for k <C 1, we are unable to predict consequences of (5.2) without 
additional explicit calculations. However, because of the late time behavior of Re^ ^> 1, we do 
not expect that this case is particularly relevant for the simulation of heavy ion collisions (see the 
more detailed argument in the paragraph below). 

In the opposite case Re^ 3> 1, Re^ <C 1, and for /t C 1, all terms containing derivatives with 
respect to y can be neglected and Eq. (5.2) can be seen as a set of equations describing fluid motion 
in 2 + 1 dimensions with rapidity entering only as a parameter. With respect to the transverse 
coordinates one would expect a turbulent flow pattern. We remark that k as defined in Eq. (5.6) 
decreases with time r and the ratio of the two Reynolds numbers (which is independent of k) 
increases with time 

Re^/Re^ = " rA y 2 . (5.10) 

Vy I 

Therefore, in contrast to the case Re^ 1, Re^ ^> 1 discussed above, the case Re^ ^> 1, 
Re {y) < 1 can persist at late times, and it can be reached dynamically. For completeness, we 
mention finally the region k ^> 1 when derivatives with respect to y enter in the damping term. 
Considered as a function of the transverse coordinates x± and x 2 , the velocity field will be damped 
by locally varying rates due to irregularities with respect to the rapidity argument. However, this 
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damping is dominated by the first non-linear term in Eq. (5.2) so that the fluid will again behave 
turbulent with respect to the transverse coordinates. 

The above discussion indicates that the case Re^ 3> 1, Re^ C 1, k <C 1 may be of particular 
relevance for the discussion of the onset of non-linear turbulent behavior in heavy ion collisions, 
since it can be reached dynamically and since it persists at late times. For a large set of initial 
conditions, we therefore expect that the fluid dynamics of heavy ion collisions evolves towards a 
system with effectively two-dimensional turbulent behavior. To address the issue to what extent the 
evolution towards turbulence could be completed within the finite duration of a heavy ion collision, 
let us finally put some numbers to this parametric discussion. Because of the late time limit of 
eq. (5.10), we focus on the transverse Reynolds number 

Re^ = V -^. (5.11) 
V 

We consider typical values for a heavy ion collision at LHC energies, e.g. T « 0.3 GeV and a length 
scale I ~ 5fm. Taking for the transverse velocity a fraction of the velocity of light, say vt = 0.1 c, 
one finds vtIT m 1 and 

Re< T > = ^-. (5.12) 
r)/s 

For small values of the normalized viscosity rj/s < 1, one therefore expects a transverse Reynolds 
number Re^ that is larger than unity but not many orders of magnitude larger than unity (Re^ < 
100). Such values are not sufficiently large to expect fully developed turbulence. A value Re^ > 1 
indicates, however, that the system can be driven outside the region of validity of a laminar evolution 
and that it may display features indicative of the onset of turbulent behavior. 



6 Qualitative features of turbulence 

In the previous section, we have discussed the general conditions under which the time-evolution 
of fluctuations in a Bjorken expansion scenario can be expected to lead to the onset of turbulent 
behavior. Despite the relativistic nature of the system under study, we found that upon coordinate 
transformation, the time evolution of fluctuations is governed by an equation that takes the form of a 
two-dimensional Navier-Stokes equation at late times. Although the Navier-Stokes equation presents 
still many deep problems, much is known about fully developed turbulence in this system. Classical 
achievements include in particular the scaling theory by Kolmogorov [59] for three-dimensional 
turbulence and its extension to the two-dimensional case mainly by Kraichnan [60] and Batchelor 
[61]. For reviews of this field, see [62] 

To the best of our knowledge, the Bjorken scenario considered here is the first example of a 
relativistically expanding three-dimensional scenario that under suitable initial conditions evolves 
dynamically into a system with effectively two-dimensional turbulent dynamics described by a non- 
relativistic Navier-Stokes equation. Turbulence in the two-dimensional case is known to display 
some characteristic qualitative differences in comparison to the three-dimensional case. Although 
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the present section contains, strictly speaking, no novel results, our finding of an effectively two- 
dimensional fluid dynamic propagation of fluctuations at late times prompts us to discuss here the 
pertinent features of fully developed turbulence that apply to the Bjorken scenario for sufficiently 
large values of Re^* 1 . In particular, we point to the phenomenon of an inverse cascade that exists 
only in the case of two-dimensional turbulence and that provides a unique mechanism for enhancing 
fluctuations on large spatial scales during the evolution towards turbulence. To set the stage for this 
discussion, we introduce first shortly a statistical description of fluctuations in the fluid velocities 
and energy densities. 

6.1 Fluctuation spectra 

In general, fluctuations in the fluid velocities and the energy densities can be described statistically 
in terms of a r-dependent probability distribution 

Pt[u^(t, xi,x 2 , y), e(r, x u x 2 , y)\ . (6.1) 

Eq. (6.1) describes an ensemble of events with equal "macroscopic" properties such as nucleon 
number, center of mass energy and impact parameter. Bjorkens model of a unique fluid velocity 
and energy density is then recovered by taking the probability distribution in (6.1) to be infinitely 
narrow. 

We consider a generalization of Bjorkens model where the assumed symmetries (boost invariance 
in the longitudinal direction and translation and rotation invariance in the transverse plane) are 
broken by the fluctuations for a particular event but hold in a statistical sense. This means that the 
probability distribution (6.1) is invariant under these symmetries. This implies in particular that 
the expectation value of velocity is given by 

(«"> = (1,0,0,0) (6.2) 

for all times r > r and for all values of xi, x 2 , y. Similarly the expectation values of thermodynamic 
scalars such as e, p, s, T etc. will be a function of r, only. In general, however, the r-dependence 
of these expectation values will differ from the ones obtained by Bjorken due to non-linear effects 
of fluctuations. 

Beyond the expectation values, the probability distribution (6.1) can be characterized in terms 
of correlation functions. In particular, the two-point correlation function of velocities at equal time 
r is defined by = 1, 2, y) 

G u( T i x i -x' v X2 -x' 2 ,y-y') = (u l (r,x 1 ,x 2 ,y) u>{T : x' 1: x' 2 ,y')). (6.3) 

Due to translational symmetries this depends only on the differences in the spatial coordinates. 
Similarly, we define 

G t (t,x 1 - x[,x 2 -x' 2: y-y') = (T(r,x 1 ,x 2 ,y) T(r, x[, x' 2 , y')) (6.4) 
and the cross-correlation 

g It( t i x i ~ x 'n x 2 ~ x 2,y~ y') = (u 3 (r,x 1 ,x 2 ,y) T(t, x[, x' 2 , y')). (6.5) 
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The generalization to other scalar quantities such as energy density e or pressure p or to un-equal 
time arguments is obvious. It is useful to introduce also the Fourier decomposition 

/d 3 k 

and the abbreviation 

G$(t) = G i i( T ,x 1 = 0,x 2 = 0,y = 0) . (6.7) 

This generalizes trivially to the other functions. We note that for symmetry reasons, one has 
G%(t) = for i ^ j and G^(t) = G 22 (t). 

For a solenoidal fluid with djVP = one has 

d i G%(T,x 1 ,x 2 ,y) = Q (6.8) 

or in momentum space 

kiG»(r,k) = 0. (6.9) 

In this framework, characterizing the fluid evolution of a heavy ion collision amounts to make 
statements about the form of correlation functions such as G l j(r, k) etc. at some given time r. It is 
clear that the form of these correlations depends strongly on the initial conditions. For example, if 
the initial fluctuations at time tq are small enough (or, equivalently, the viscosities large enough to 
have small Reynolds numbers) so that the linearized equations (4.3) - (4.6) can be applied, one can 
derive from them linear evolution equations for the set of correlation functions G%, Gj and G J ud . 
The form of these functions at time r is then directly linked to the corresponding functions at time 

TO- 

The situation is much more complicated in the presence of non-linear contributions to time 
evolution, even if the system is far from the conditions of fully developed turbulence. Indeed, if 
one tries to use the non-linear equations (3.14) and (3.15) to derive evolution equations for the 
two-point correlation functions, one finds that three-point correlations get involved, as well. The 
evolution equation for these involve even higher correlations and so on. This is an instant of the 
well-known closure problem in the statistical description of fluids. The same problem appears in non- 
perturbative formulations of quantum and statistical field theories. No exact analytical solutions 
are known and advanced techniques from statistical mechanics and field theory are needed to find 
approximate ones. The scaling theory of Kolmogorov provides some insight into these problems. 

6.2 Scaling theory of turbulence at large Reynolds number 

Any fluctuation present in the initial conditions of a relativistic heavy ion collision can be regarded 
as a source of non-thermal, say 'mechanical' energy. Most generally, one would like to understand 
to what extent this mechanical energy dissipates to thermal energy, and to what extent it does 
not but leaves characteristic, dynamically evolved fluctuations visible in final state observables. 
In section 4, we have shown examples of the dissipation (or amplification) of some fluctuation 
modes in a linearized description that applies to very small Reynolds numbers. Here, we want to 
comment on the opposite case of very large Reynolds number or very small viscosity. For large 
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Reynolds numbers, we have found in section 5 that fluctuations on a Bjorken background field 
are governed by an evolution equation that takes for sufficiently late times the form of a non- 
relativistic Navier-Stokes equation. Also, for situations of sufficiently small Mach number, there 
is a rationale for setting 9 — > 0. Therefore, our discussion of the case of large Reynolds number 
will reduce to recalling pertinent features of Kolomogorov's theory of homogeneous turbulence for 
a non-relativistic, incompressible fluid at very large Reynolds numbers. 

For a non-relativistic fluid, it is common to denote the fluid kinetic energy per unit mass by 
^[v 2 ], where the square brackets stand for spatial averaging. The rate of dissipation to thermal 
energy for a three-dimensional evolution is given by 

d X 1 

£diss = -^2^ = 2 U ^^( diVj + djVi ^ = u ^ ' ( 3 - dim - case ) (6-10) 

This shows that energy dissipation is mainly due to fine structures of the velocity field for which the 
gradients are large; dissipation is mainly taking place at large wave-vectors k. If one now decreases 
the viscosity, one finds finer and finer structures emerge so that the energy dissipation rate Edi SS 
remains positive. In fact, the mean-square vorticity \[ui 2 ] (also called enstrophy) grows ~ 1/v for 
v — y 0. This is possible since \[ti 2 ] is not only given by the vorticity present at some initial time or 
generated from an external driving force. Rather, it can be generated also by non-linear terms in 
the three-dimensional Navier-Stokes equation. The mechanical energy is cascaded from the large 
length structures to the smaller ones by virtue of the non-linear terms in the Navier-Stokes equation. 
This is the famous cascade picture of Richardson [63]. In his words: "Big whorls have little whorls, 
Which feed on their velocity; And little whorls have lesser whorls, And so on to viscosity." 

Let us now come to the situation in two spatial dimensions. The most important difference 
to the three-dimensional case concerns the evolution of mean-square vorticity in the absence of an 
external driving force (vorticity u = d\V2 — c^Ui has now only a single component) 

~[uj 2 } = -v\(yu) 2 }. (2-dim. case) (6.11) 

This shows that \[oj 2 } never increases as a function of time. This in turn implies that energy per 
unit mass is conserved for vanishing vorticity, 

~\v 2 ]^0 for 
dtT J 

A cascade of mechanical energy from large structures into smaller ones where it is finally dissipated 
is therefore not possible 6 . 

Without going further into the detailed mechanism of two-dimensional turbulence let us now 
attempt to transfer some of the insights that have been gained in this field to heavy ion physics, in 
particular fluctuations around Bjorken flow. We consider the case of large fluctuations and small 

6 It has been argued that a cascade can take place for enstrophy \\^ 2 \ instead of energy ^[u 2 ], however. This is 
possible if [(Vw) 2 ] grows ~ 1/v for v — > due to non-linear terms. 
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kinematic viscosity v so that the Reynolds number Re is large. Also, we concentrate on the limit of 
large time r where (5.2) becomes two-dimensional. Similar to the non-relativistic case one has now 

j[vl + vl\^U for u^O. (6.13) 

For a fixed time t and rapidity y one can characterize the correlations of the velocities in the 
transverse plane by (m, n — 1, 2) 

(G v ) mn (t, xi — x 1 , x 2 — x 2 , 0) = (v m (t,xi,x 2 ,y)v m (t,x' 1 ,x' 2 ,y)) 

-m\2/3 (6.14) 
G™ n (T,t,x 1 -x , 1 ,x 2 -x' 2 ,0). 

Adapting to the standard notation used in the literature about turbulence, we write the Fourier 
transform of this as 



X!,x 2 ,y) 



r d 2 k i(klX1+k2X2) k m k n \ 2n E(f k) (6 m) 

J (2vr)2 e kl + kl) k h[t ' k) - {bAb) 

The tensor structure of (6.15) follows from rotational invariance and from the solenoidal constraint 

(5.3). The function E(t,k) depends on k\ and k 2 only in the combination k = ykf + k 2 . The 
normalization in (6.15) is chosen such that 

1 1 2 r°° 

E(t) = -(v 2 1 +v 2 2 ) = -J2(Gv) mm (t,0,0,0)= / dkE{t,k). (6.16) 

™ i J 



For a non-relativistic fluid, the function E(t, k) describes how the fluid kinetic energy per unit mass 
is distributed over the different wave vectors. We emphasize that in the relativistic setup considered 
here, E(t, k) is not directly representing kinetic energy. Instead it simply parameterizes the con- 
tribution to the fluctuating transverse velocity field from different wave vectors. The contribution 
of these fluctuations to kinetic energy, for example in the laboratory frame, can be determined but 
the resulting relation is more-complicated than in the non-relativistic case. 

For the case of a relativistic heavy ion collision, the initial distribution E(t , k) would charac- 
terize the relative strength with which different length scales 1/k are represented in the fluctuating 
initial conditions. The question of how this distribution of fluctuations evolves amounts then to 
studying the time-dependence of E(t,k). Here, we point only to one remarkable feature of the 
time-dependence of a two-dimensional fluid at large Reynolds number, that can be understood in 
the scaling theory of freely decaying turbulence in two dimensions developed by Batchelor [61]. 
This theory is based on the assumption that at sufficiently late times, the function E(t, k) remem- 
bers only a single number from its initialization, namely its average fluid velocity A defined by 
A 2 = \ {v\ + v 2 ). From dimensional reasoning it follows then that 

E(t,k) = XHh(kXt). (6.17) 

It follows from (6.16) that the dimensionless function h(x)is normalized to unity, J °° dx h(x) = 1. 
Interestingly, if one assumes that E(t, k) is dominated by the region around some characteristic 
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wave vector ku-, then this scale will change with time according to 




(6.18) 



This implies that kinetic energy is shifted from small length scales to larger ones, in contrast to the 
Richardson cascade in three dimensions [63]. We would like to close this section on a cautious but 
speculative note: We recall first that Batchelor's theory of freely decaying turbulence was developed 
for very large Reynolds numbers that may not be realized in heavy ion collisions. However, the above 
considerations may make it conceivable that non-linear effects in the fluid dynamic evolution related 
to the onset of turbulence can move fluctuations in the initial kinetic energy to larger spatial scales. 
This phenomenon would be a distinct characteristics of an effectively two-dimensional turbulent 
evolution, and it may be identified experimentally by finding fluctuations related to length scales 
that are inconceivable to be present in fluctuating initial conditions. 

7 The sensitivity of particle spectra on velocity correlation functions 

In the previous section 6, we have seen that the fluid dynamic evolution of fluctuations can be 
characterized efficiently in terms of correlation functions of fluid velocities. Here we discuss how 
information about such velocity correlations enters the particle spectra that are experimentally 
accessible in heavy ion collisions. 

The starting point of our discussion is the 'freeze-out' phase space distribution f(x,p) that 
parametrizes the matter distribution at the time Tf when the particles decouple from the fluid 
dynamic evolution. In the following, we shall view f(x,p) as describing an event-averaged smooth 
fluid system supplemented by event-specific fluctuations. We shall then ask how these fluctuations 
are reflected in observables. This logic should apply to arbitrary choices of f(x,p). For the purpose 
of illustration, however, we shall restrict our discussion to a simple ansatz 7 that describes a locally 
approximately thermal distribution consistent with the Bjorken background field of section 3, and 
that allows for the implementation of local fluctuations on top of this background field. A simple 
choice with these properties is the Boltzmann distribution 



where the normalization d is fixed by the spin and flavor degeneracy of the degrees of freedom that 
decouple from the system. Assuming that the freeze-out takes place at some proper time Tf Q when 
the average temperature drops below some freeze-out temperature T fo , the hadronic spectra can be 
calculated using the Cooper-Frye freeze-out prescription 



7 We note as an aside that phenomenologically more realistic choices of f(x,p) would be significantly more com- 
plex. In particular, they would contain information about the finite spatial extent of the matter distribution in the 
transverse and longitudinal distribution, they would supplement the ideal gas expression (7.1) by terms proportional 
to viscosity [12-14], and they may implement correct Fermi-Dirac or Bose-Einstein statistics instead of the Boltzmann 
distribution (7.1). 



p M it M (x) 



f(x,p) = de n*) 



(7.1) 




(7.2) 
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Here, the freeze-out volume is determined by p^dY,^ = mTTcosh(r] — y)dx\dx 2 dy for the case of 
Bjorken expansion. 

In practice, the spectra (7.2) measured in heavy ion collisions include averaging over many 
events. On the level of the freeze-out distribution f(x,p), this event average corresponds to an 
ensemble average with respect to the fluid velocity, energy density or temperature fields, respectively. 
Denoting the corresponding averaging by triangular brackets, we replace therefore in the calculation 
of (7.2) the function / by 

f(x,p)=d(e V (7.3) 

Denoting the particle four- momentum by (p° , p 1 , p 2 , p 3 ) = (jht cosh y , p 1 , p 2 , ttit sinh y) with trans- 
verse mass squared m 2 - = p\ + m 2 , p\ = (p 1 ) 2 + (p 2 ) 2 , we expand f(x,p) for small fluctuations 
around the velocity profile of Bjorken (u°, it 1 , u 2 , u y ) = (cosh rj, 0, 0, sinh rj) and the constant freeze- 
out temperature Tf 

rriT cosh(r/ — y) 



f(x,p) = de 



vrirp coshf^ — y) 
Wo 



1 + 



T 2 

J fo 



-{T - T fo ) + — (pxu 1 + p 2 u 2 + rm T sinh(?7 - y) u y ) 

J- to 



m\ cosh 2 ^ — y) itlt cosh(?y — y) 



+ ^2 (piu 1 + p 2 u 2 + t m T sinh(r] - y) u v ) 2 + ( 

^fo V 

H t^a ipiu +P2U + Tm T smh[T] - y) u y )(T - T fo ) + . . . ). 
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The terms linear in T — Tf Q or vanish by definition or due to symmetry reasons and similar the 
cross-terms ~ u^u 1 with % ^ j. Also, due to translational symmetry the variances at r = Tf are 
actually independent of the coordinates x\,x 2 and y, such that 

((u y ) 2 ) = Gf (r fo ), (7.5) 
((T-T io ) 2 ) = G T (r io ). 

This allows one to write the one-particle spectrum for small fluctuations around Bjorken flow as 
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In the simple model studied here, information about velocity correlations enters the spectrum only 
via the coordinate-independent three numbers G^(t{ ), G^(t{ ) and Gt(t{ ). For another model 
choice, the information may be slightly different. For instance, if one would replace the sharp 
freeze-out at r fo by a decoupling at times r around r fo , then the three numbers G^(t {o ), G^ y (r {o ) 
and Gt{t{ q ) would be replaced by averages over r. Irrespective of such model-dependent nuances, 
however, it is a generally known feature that the single particle spectrum (7.2) is only sensitive to 
space-time averages over the distribution f(x,p) and therefore does not contain information about 
correlations between different space time points. 

7.1 Generalization to identical two-particle correlations 

As discussed in section 6, the dependence of velocity correlations on the wave- numbers ki,k 2 and 
k y allows for a detailed characterization of fluid dynamic behavior, including information about the 
dissipation of fluctuations and the manifestations of turbulence. The one-particle spectra discussed 
so far contain only the information (7.5) about the correlations of fluid fields 

G*(T,x x ,x 2 ,y), G 3 uT (T,x u x 2 ,y), G T (r,x 1 ,x 2 ,y) (7.8) 

at equal positions (xi = x 2 = y = 0) . Here, we point out that identical (Bose-Einstein) two-particle 
correlation functions are linear functionals of (7.8) and may thus provide information about the 
wave number dependence of velocity correlations. 

Two-particle spectra for pairs of identical bosons (sb/f — 1) or fermions (sb/f — — 1) of 
4- momenta Pa, Pb respectively, can be written as [64] 



E A E B „ dN „ — = f (pAl^dE!* {pB) v dT! v f(x,p A )f(x',p B ) 
d 6 p A d 6 p B J 

+ SB/F J \{p A +P B ),dW I(p A +p B )^ e *(PA-PB)M(*-^ / ( Xj E^) / ( x / ) £^). 



(7.9) 



Here, the first term is what one would expect from kinetic theory for classical particles while the 
second term results from the quantum statistics of identical particles. Data about two-particle spec- 
tra are typically normalized by a mixed-event technique that corresponds to forming a normalized 
correlation function 



C(pa,Pb)= "Z^ - (7.10) 
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As it stands, equation (7.9) is valid for an event-specific realization of the velocity and tempera- 
ture fields. Experimental data are for event samples that correspond to averages over the hydrody- 
namic fields. This amounts to replacing in equation (7.9) the products of phase-space distributions 
/ by the corresponding event averages {f (x ,p a) f (x' ,pb)) an d (f(x, PA 1 2 PB )f(x', PA 1 2 PB )). Paralleling 
the arguments employed for the calculation of the one-particle spectra via (7.4), one should then 
expand the arguments of {f (x, pa) f (x' , Pb)) for small fluctuations around the event-averaged back- 
ground fields. In general, the arguments of these averages depend on space-time difference x — x', 
and they contain information about the relative position dependence of the correlations (7.8) in the 
fluid fields. 
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In practice, the way in which information about (7.8) enters the two-particle correlation func- 
tions may depend significantly on model-specific choices for f(x,p). For illustrative purposes, we 
explore here the particularly simple Bj or ken-like model without spatial constraints in the transverse 
direction. We ignore all correlations involving the temperature field on the ground that these are for 
a compressionless situation formally of higher order in the fluctuating velocities, see Eq. (5.4). For 
the discussion in the following, we also neglect the rapidity-dependence of the correlation functions 
thus eliminating many terms proportional to G % £ with i = y or j = y (or both), as well as terms 
involving G v uT (The assumption of a vanishing rapidity dependence is relaxed in Appendix A). With 
these approximations, we concentrate therefore on the correlation functions of velocities G™ n with 
m, n = 1, 2. Due to rotational invariance in the transverse plane, the velocity correlation in Fourier 
space can be written in the form 

ran 9lij~i k m k n 92(r,k)] , (7.11) 

with k = \fk\ + k\. In terms of the pair momentum P = \(pa + Pb) and the relative momentum 
q = pa ~Pb, the two-particle correlation function at mid-rapidity T]a = t]b = takes then the form 
(see appendix A for details of the derivation) 

C(P,q) = 1 + S B/F {2v)W\q T )±- + P| rjf /4 <7l(T fo ,0) 
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(7.12) 



with P T = (Pi, P 2 ), q T = (qi, qi) and q T = \fq%. 

Realistic phase space distributions f(x,p) have support in a finite transverse area At only, 
and this would lead to a fall-off of the correlation function to unity in the relative transverse 
momentum ^ on a scale of order 1 / y/ At- Often, this fall-off is parametrized by a Gaussian ansatz 
in terms of HBT radius parameters, so that the first two terms of (7.12) would take the form 
C(P,q) = 1 + sb/f exp [—At q T }- For the simplified model of infinite transverse extension discussed 
here, this contribution is singular oc (27i) 2 S^ 2 \qr) and we have written it in a formal way normalized 
to unity at q\ = 0. 

For the purpose of the following discussion, the transverse translational invariance of the present 
toy model presents the technical advantage that the gr-dependence of the correlation function 
(7.12) is solely dependent on the wave number dependence of velocity correlations. Effects that 
could confound the interpretation of the (f^-dependence in practice, such as effects from a finite 
transverse geometry and from transverse (event-average) velocity gradients, are not included in 
the present model. Therefore, the following discussion allows us to illustrate how information 
about velocity correlations enters two-particle correlation functions, but it limits our discussion of 
how such information could be disentangled from other effects in a phenomenologically relevant 
scenario. Keeping this caveat in mind, we observe that the first term in (7.9) can be viewed as an 
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incoherent superposition of single-particle spectra and therefore does not furnish information that 
is not yet contained in single-particle spectra. In contrast, the quantum-statistical second term 
oc sb/f furnishes novel information about the wave number dependence of <7i(rf , q?) and ^2(^0, Qt)- 
We note as a curious aside that for a situation of fully developed two-dimensional turbulence, 
one can predict the form of gi(r, k) and ^(t, k) from the scaling theory of Kraichnan and Batchelor. 
In particular, comparing (6.14), (6.15) and (7.11) one can express g% and gi as a function of E(t, k). 
Moreover, from Batchelors scaling theory of freely decaying two-dimensional turbulence one finds 
then the following scaling in the inertial range 

c 



gi{T,k) - _ |n H;J . 

(7.13) 



r 10/3 k 4 

02 (r, k) = 



r 10/3£.6 ' 

with a common but unpredicted constant c that reflects the absolute scale of velocity fluctuations. 
For the correlations function (7.12), this results in an additive term with a very slow power-law 
fall-off of the form oc A ^ j, 2 g4 . Interestingly, if the velocity correlations occurs on scales that are 

significantly smaller than the transverse extension a/ At of the system, then the slow power-law q?- 
dependence persists at relative momentum scales that are significantly larger than the typical scales 
1/R set by HBT radius parameters. It is an exciting possibility that a measurement of a power-law 
l/g^.-dependence in two-particle correlation functions may provide a characteristic signature for 
turbulent distributions in the fluid dynamically evolved velocity fields of a heavy ion collision. We 
caution that the model discussed here is a simplified one; also, for intermediate Reynolds numbers 
one expects corrections to the case of fully developed turbulence, see e.g. [62]. What may persist 
in a phenomenologically realistic scenario, however, is the general idea that velocity correlations on 
small scales It induce two-particle correlations on large scales qT ~ and that these correlations 
are expected to be governed by a power-law fall-off. 



8 Discussion and Conclusion 

Recent data analyses from RHIC and LHC have given support to arguments that soft hadron spectra 
may result from the fluid dynamic response to initial conditions with significant event-by-event 
fluctuations. Motivated by this suggestion, we have studied here how event-by-event fluctuations 
propagate on top of an event- averaged fluid dynamical background of Bjorken type. The choice 
of a Bjorken background field is a simplification that retains essential elements of the expected 
fluid dynamical evolution of relativistic heavy ion collisions. As shown in the present paper, it 
allows for a particularly explicit, partly analytical discussion of the propagation of fluctuations in 
an expanding fluid dynamic system. In particular, we have found for the case of a laminar evolution 
explicit expressions for the attenuation or amplification of all fluid dynamic modes over the time 
scale relevant in heavy ion collisions. And we have specified the general conditions for non-linear 
effects in the dynamics of fluctuations, finding in particular that the late time dynamics evolves 
towards an essentially two-dimensional system with turbulent behavior. Here we discuss our main 
findings in more detail: 
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To discuss the propagation of fluctuations, one needs to specify first the nature of the fluc- 
tuations that are propagated. Recent studies have focussed mainly on fluctuations in the energy 
density (or, equivalently, entropy density), where the Glauber model provides a phenomenologically 
supported basis for assuming local fluctuations on a particular transverse scale. However, it had 
been pointed out already that fluctuations in the velocity field may be present in the initial condi- 
tions for fluid dynamic evolution. Velocity fluctuations could arise from pre-equilibrium evolution, 
or they could be a natural consequence of fluctuations in the primary interactions of elementary 
constituents. To the best of our knowledge, there is no a priori argument that fluctuations in energy 
density dominate over fluctuations in other fluid dynamic fields. Also, it requires studies allowing 
for all possible fluctuations to address the question whether and how fluctuations in velocity and 
energy density can be disentangled. In section 2, we argued that a discussion of the fluid dynamic 
response to fluctuating initial conditions should be based on a formulation that allows for fluctu- 
ations in all fluid dynamic fields. In particular, we supported with a model study the idea that 
fluctuations in the velocity fields may carry significant vorticity. For a fluid with conserved charges 
such as baryon number and electric charge one should take fluctuations in these quantities into 
account, as well. 

In general, a separation of fluid dynamical fields into background and fluctuations is not neces- 
sary. For instance, recent studies of event-by-event fluid dynamics propagate event samples of initial 
conditions numerically without separating fluctuations from background fields. In principle, such 
full fluid dynamical simulations allow to explore under the most versatile model assumptions the 
fluid dynamic response to fluctuating initial conditions. But an explicit mode-by-mode formulation 
of the dynamics of fluctuations around a fluid background field seems well-suited to study which 
fluctuating modes in the initial conditions can survive the strong dynamical evolution in a heavy 
ion collision unattenuated, whether there are mechanisms that may amplify some modes, and which 
modes are 'filtered out' by the medium due to dissipative effects in the fluid dynamic evolution. To 
address such questions in a simplified framework that accounts for the main features of fluid ex- 
pansion, we have formulated in section 3 the fluid dynamical evolution of local fluctuations around 
average fluid fields of Bjorken type. For this system, we have found a peculiar rescaling of the time 
variable, t oc r 4//3 , that allowed us to write the relativistic fluid dynamic evolution of fluctuations 
in a form resembling a non-relativistic Navier-Stokes equation. 

If the Reynolds number of a fluid system is not too large, then important elements of the 
dynamics may be understood in a linearized treatment. In section 4, we observed that in this laminar 
case, fluctuations in vorticity decouple from sound modes and fluctuations in energy density. For a 
parameter set that is characteristic for heavy ion collisions, we have then studied how different modes 
are amplified or attenuated over the time scale of order 10 fm/c of a heavy ion collision. Remarkably, 
the longitudinal vorticity modes have an algebraic enhancement factor that can overcompensate on 
this time scale the typical exponential decay due to dissipative effects. Therefore, vorticity, if 
not present in the initial conditions, will not be generated as long as the dynamical evolution 
is laminar, and it is therefore unlikely to be generated in a sizable amount for small Reynolds 
numbers. However, if present in the initial conditions, some vorticity modes will be amplified 
significantly during the dynamical evolution. To the best of our knowledge, this phenomenon has 
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not been studied yet in numerical simulations, and it would be interesting to see how it manifests 
itself in the presence of other background fields. In addition to the vorticity modes, we have also 
explored the propagation of sound modes that result from initial fluctuations. In general, we find 
that fluctuations of sufficiently long wave-length pass unattenuated over time scales relevant for 
heavy ion collisions, while short wave-lengths that reflect finer structures in the fluctuating initial 
conditions, are dissipated on shorter length scales. There is also a characteristic difference between 
transverse components, and the components that propagate in the longitudinal direction in which 
the system expands according to Bjorken's model. 

Outside the regime of validity of a linearized treatment, the discussion of solutions of fluid 
dynamics is very complicated and typically requires numerical techniques. For our problem of 
fluctuations around Bjorken flow, we are in the special and fortunate case that we can relate 
the full relativistic dynamics of fluctuations in rescaled coordinates to a non-relativistic Navier- 
Stokes equation. This allows us to discuss the possibility of a turbulent evolution in terms of those 
concepts and parametric estimates that have proven useful in characterizing turbulent phenomena 
of non-relativistic systems. In section 5, we introduce both a longitudinal and a transverse Reynolds 
number to characterize the non-linear dynamics of fluctuations on top of a Bjorken background field 
that shows strong dynamical expansion only in the longitudinal direction. We find in particular 
that the late time dynamics will evolve a large set of initial conditions into a regime where the 
dynamical evolution is effectively two-dimensional, and where the transverse Reynolds number can 
be sizable. This indicates a window for a two-dimensional, non-linear evolution towards turbulent 
behavior. Motivated by this observation, we have summarized in section 6 characteristic features 
of turbulent behavior, detailing the differences between the cases of two-dimensional and three- 
dimensional evolution. We recall from this discussion in particular that in three dimensions, fluid 
kinetic energy thermalizes typically by dissipating into vorticity modes of increasing wave number, 
i.e. decreasing length scales. As first observed by Kraichnan, this mechanism is not possible for a 
two-dimensional fluid system where one finds an inverse cascade: kinetic energy gets propagated to 
larger length scales. 

The transverse Reynolds numbers estimated in section 5 do not support the assumption that 
heavy ion collisions create systems with fully developed turbulence. However, our discussion in sec- 
tion 5 shows that realistic Reynolds numbers are not small enough to neglect non-linear effects in the 
dynamical evolution. Therefore, while we have no rationale to expect that the correlation functions 
of fluid dynamic fields generated in heavy ion collisions satisfy the scaling laws of Kolmogorov's 
theory of fully developed turbulence, we do expect that a non-linear dynamics that can be regarded 
as the onset of turbulent behavior may result in interesting (power-law) wave-number dependences 
of correlations of fluid fields. By supplementing a standard blast-wave model with event-by-event 
fluctuations in fluid fields, we have established in Sect. 7 how such fluid dynamic correlation func- 
tions manifest themselves in one- and identical two-particle spectra. For one-particle spectra, we 
find that fluctuations are a confounding factor in interpreting the m^-dependence of spectra in a 
fluid dynamic scenario. For two-particle correlations, we observe a dependence that may provide an 
experimentally accessible signature for the onset of turbulence in heavy ion collisions. The obser- 
vation is that identical two-particle correlations at large relative momentum are sensitive to spatial 
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scales that are much smaller than the transverse size of the particle producing source. If two-point 
velocity correlations show turbulent behavior on these small scales then this will translate into a 
characteristic power-law tail of identical two-particle correlations at large relative momentum. 

Let us close with a short outlook of how observations made in this study could be pursued 
further. Our study was largely motivated by the question of how initial conditions with significant 
velocity fluctuations (in the solenoidal an in the irrotational part) propagate fluid dynamically. 
Here, full fluid dynamic simulations including initial velocity fluctuations could provide further in- 
sight, for instance by evolving fluctuations around average fluid fields with more realistic transverse 
dependencies, and by quantifying to what extent initial velocity fluctuations could contribute to 
the observed azimuthal asymmetries in momentum space. Full fluid dynamic simulations including 
initial velocity fluctuations could also allow for a detailed characterization of how the scale depen- 
dence of fluid correlation functions of the type (7.8) builds up during the fluid dynamic evolution. 
This would provide in particular insight into the question on which time scales and over which range 
of wave vectors correlation functions of fluid fields may develop power-law dependences that can be 
regarded as precursors of turbulent phenomena 8 . Also, the semi-analytical approach used in the 
present work may be pursued further. In the present work, we have seen how single vorticity modes, 
and sound modes are amplified or filtered out by the dynamical evolution, depending on their wave- 
number. We plan to investigate, whether a similar understanding can be gained for a more realistic 
background field by expanding fluctuating fluid fields in terms of appropriate sets of functions so 
that one can calculate explicitly how the fluid dynamic evolution mixes different components in the 
evolution. We expect that such studies could provide an intuitive understanding e.g. of the time 
scales on which density fluctuations feed sound waves, or on which vorticity modes cascade to other 
scales. This may help significantly in the interpretation of the complex fluid dynamic phenomena 
that we expect to find realized in heavy ion collisions. 

A Explicit expressions for the two-particle spectrum 

Here, we provide further details about the calculation of the two-particle correlation function (7.12), 
and how the calculation of this correlation function could be generalized to include the rapidity 
dependence of velocity correlations, as well as effects of temperature fluctuations. In general, we 
represent fluctuations in the fluid dynamic fields in Fourier space according to 



We note in this context that the time scale on which non-linear contributions start to matter in the evolution of 
fluctuations may depend sensitively on the size and scale of the fluctuating initial condition. 




(A.l) 



33 



For the first term of the two-particle spectrum (7.9), we find then the following contribution (m^, 
mf, t]a, and t]b are the transverse masses and rapidities of particles A and B) 
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For an infinite extension in the transverse plane it is easy to perform the integrals over xi,x 2 ,Xi, x 2 . 
The resulting Dirac distributions can be used to perform the integrals over k±, k 2 , k[, k' 2 . Also, the 
integrals over y and y' can be done analytically. Transforming then back to Fourier space, one finds 
for the first term of (7.9) the expression 
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where the functions denote linear combinations of Bessel functions of the second kind, 

E (x, q) = K 1+iq (x) + Ki„ iq (x), 
Ei(x,q) = -K 2+iq {x) + -K 2 _ iq {x) + K iq (x), 

E 2(x,q) = ^K 2+iq (x) - -K 2 - iq {x) . 
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For k y = 0, this expression reduces to the first two terms in the first line of (7.12). 

Let us now consider the second term ~ sb/f in Eq. (7.9). This term is of the form of a single 
integral over the freeze-out volume 
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(A.5) 



times its complex conjugate. Expanding it up to terms that are quadratic in the fluctuations, 
one obtains two distinct contributions, In one case, the bilinear terms in the fluctuating fields are 
written at different points x, x', in the other case they are both taken at the same point. We 
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consider both cases separately. We work in the following with average pair momentum P M = 
|(Pa + Pb) = (fn^cosh-VAB, P 1 , P 2 , Tn^ B smh.T]AB) and we define a transverse mass defined m^ B = 
a/ (P ) 2 — (P 3 ) 2 and a rapidity t\ab = arctanh(P 3 /P°). This allows us to write q^x^ = — [m^cosh^^- 
y) - m B cosh(r] B - y)]r + q T x T with q T = (q l: q 2 ) and f T = (xi, x 2 ). 

The contribution to the second term in Eq. (7.9) that is quadratic in fluctuations at the same 
space-time point, can then be written as 

d Tf m^ B 
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Here, the integral over the transverse coordinates leads to a factor (27r) 2 £( 2 ) (&) (where (27r) 2 <5( 2 ) (0) = 
At is understood). This is a consequence of the fact that in the present model the transverse ex- 
tension of the particle emitting source is not limited. At mid-rapidity, rj A — i]b — Vab = 0, the 
term (A.6) simplifies to 
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With the approximations used in section 7, this expression reduces to the last terms in the first line 
of (7.12). 

We now turn to the contribution ~ sb/f where the fluctuating fields have different space-time 
argument. This term is of similar structure as (A. 2) and reads 
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The integrals over the transverse coordinates x\, x 2 , x[, x 2 leads to factors {2n) 2 8^ 2 \qT + hr) and 
(27i) 2 S ( - 2 \qr — k' T ), respectively. These can be used to perform the integrals over the transverse 
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components of k and k'. For simplicity we concentrate again on particles at mid- rapidity, t]a = 
Vb = Vab = 0. It is then straight-forward to perform the integrals over y and y' In terms of the 
fluid dynamic correlation functions, the result can then be written as 
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(A.9) 



We note that (A.9) contains information about the full fluid dynamic correlation functions in mo- 
mentum space. Under the assumptions made in section 7, this expression reduces to the last line 
of (7.12). 
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